This notebook contains the code of the paper "Bayesian Calibration of Imperfect Computer Models using Physics-Informed Priors". The models are fitted in rstan and the code is available in the folder "STAN/Approximations".

Load packages

knitr::opts_chunk$set(
  message=FALSE,
  warning = FALSE,
  comment = '', 
  fig.width = 6, 
  fig.height = 4,
  fig.align = 'center'
)
getwd()
[1] "/Users/michais/Desktop/BCPI_codes_new"

Reality and modelling choice

\[\begin{align} \mathcal{R}: \quad & \frac{d P(t)}{d t} + \frac{P(t)}{R_2C} = \frac{Q(t)}{C} \left (1 + \frac{ R_1}{R_2} \right ) + R_1 \frac{d Q(t)}{dt} \quad \text{ (the misspesified model we use to fit the data) }\text{ [WK3] }\\ \eta: \quad & Q(t) = \frac{1}{R}P(t) + C \frac{dP(t)}{dt} \quad \text{ (the model we use to simulate data) } \text{ [WK2] } \end{align}\]
# uncomment to install
# install.packages("rstan")
# install.packages("ggplot2")
# install.packages("tidyverse")
library(rstan)
library(ggplot2)
library(tidyverse)
rstan_options(auto_write = TRUE)
options(mc.cores = 3) # allocate 3 cores (for each model we run 3 chains in parallel)
# Numerical simulator of the WK3 model
source("functions/WK2and3_sim_fn.R")
# Load flow data 
d = readRDS("Data/Inflow_time.rds")
Rtrue = 1; Ctrue = 1.1; Ztrue = 0.05 
flow = d$inflow*0.95
time = d$time
nP = 90 # number of pressure data
nI = 100# number of inflow data
nc = 1  # number of cardiac cycles
nflow = length(flow)
# 1. simulate WK3 data (R=R_2, Z=R_1)
Psim = WK3_simulate(flow = flow, time = time, R = Rtrue, C = Ctrue, Z=Ztrue) # simulate WK3 data for a given flow Q(t)
P_true = Psim
# 2. choose pressure and inflow indices
indP = round(seq(1, nflow, length.out = nP)); indI = round(seq(1, nflow, length.out = nI))
yP_real = Psim[indP]; yI_real = flow[indI] # noise free fimulated pressure and flow
# 3. Add noise
# set.seed(0)
set.seed(1)
Pnoise = rnorm(nP*nc, 0, 4) # sample pressure noise from N(0, 4^2)
Inoise = rnorm(nI*nc, 0, 10) # sample flow noise from N(0,10^2)
yP_real = rep(yP_real,nc) 
yI_real = rep(yI_real,nc)
# 4. store individual data in the population matrices
yP = yP_real + Pnoise # add noise
yI = yI_real + Inoise # add noise
tP = time[indP] # corresponding time (synchronized for the two cycles)
tI = time[indI] # corresponding time (synchronized for the two cycles)
data_PI = list(nP=nc*nP, nI=nc*nI, tP=rep(tP,nc), tI=rep(tI,nc), yP=yP, yI=yI, mP=12, mI=10)
WK2_VFE = stan_model('STAN/Approximations/VFE/WK2_delta_VFE.stan')
kp = kmeans(data.frame(x=data_PI$tP), centers = data_PI$mP)
ki = kmeans(data.frame(x=data_PI$tI), centers = data_PI$mI)
init = list("zP" = as.vector(kp$centers), "zI" = as.vector(ki$centers))
op_VFE=optimizing(WK2_VFE, data=data_PI, hessian=FALSE, init = init, verbose=TRUE, seed=0)
Chain 1: Initial log joint probability = -13400
Chain 1:     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
Chain 1:       19      -962.859     0.0138579        2.3978           1           1       31   
Chain 1:     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
Chain 1:       39      -949.581     0.0272867       2.43968           1           1       56   
Chain 1:     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
Chain 1:       49      -949.485   0.000188747     0.0088334           1           1       68   
Chain 1: Optimization terminated normally: 
Chain 1:   Convergence detected: relative gradient magnitude is below tolerance
op_VFE
$par
         rho        alpha        rho_d      alpha_d       mu_wk2       sigmaP       sigmaI            R            C        zP[1]        zP[2]        zP[3]        zP[4] 
 0.195718978  4.897341220  0.057918221  0.001593782 97.398040198  9.760864357 63.624564498  0.955552758  1.013489858  0.080072408  0.133867715  0.455149871  0.194511717 
       zP[5]        zP[6]        zP[7]        zP[8]        zP[9]       zP[10]       zP[11]       zP[12]        zI[1]        zI[2]        zI[3]        zI[4]        zI[5] 
 0.274922964  0.724878480  0.826984618  0.944283254  0.365017283  0.544772397  0.020043511  0.634880324  0.970363221  0.614334337  0.714632658  0.049797442  0.394214611 
       zI[6]        zI[7]        zI[8]        zI[9]       zI[10] 
 0.511302773  0.809613912  0.274381891  0.161115333  0.894987597 

$value
[1] -949.4846

$return_code
[1] 0

$theta_tilde
          rho    alpha      rho_d     alpha_d   mu_wk2   sigmaP   sigmaI         R       C      zP[1]     zP[2]     zP[3]     zP[4]    zP[5]     zP[6]     zP[7]     zP[8]
[1,] 0.195719 4.897341 0.05791822 0.001593782 97.39804 9.760864 63.62456 0.9555528 1.01349 0.08007241 0.1338677 0.4551499 0.1945117 0.274923 0.7248785 0.8269846 0.9442833
         zP[9]    zP[10]     zP[11]    zP[12]     zI[1]     zI[2]     zI[3]      zI[4]     zI[5]     zI[6]     zI[7]     zI[8]     zI[9]    zI[10]
[1,] 0.3650173 0.5447724 0.02004351 0.6348803 0.9703632 0.6143343 0.7146327 0.04979744 0.3942146 0.5113028 0.8096139 0.2743819 0.1611153 0.8949876

FITC

zP_opt_VFE=op_VFE$par[grep("zP",names(op_VFE$par))]
zI_opt_VFE=op_VFE$par[grep("zI",names(op_VFE$par))]
# plot(sort(zP_opt_VFE))
# plot(sort(zI_opt_VFE))
data_PI_Z_VFE= data_PI
data_PI_Z_VFE$zP = zP_opt_VFE
data_PI_Z_VFE$zI = zI_opt_VFE
fit_post_VFE=stan(file='STAN/Approximations/VFE/WK2_delta_VFE_fixed_Z.stan',
                  data=data_PI_Z_VFE,
                  chains=3,
                  iter=1000,
                  seed=0
)
starting worker pid=67458 on localhost:11354 at 11:30:29.975
starting worker pid=67472 on localhost:11354 at 11:30:30.172
starting worker pid=67486 on localhost:11354 at 11:30:30.365

SAMPLING FOR MODEL 'WK2_delta_VFE_fixed_Z' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.009392 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 93.92 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'WK2_delta_VFE_fixed_Z' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.010761 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 107.61 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 

SAMPLING FOR MODEL 'WK2_delta_VFE_fixed_Z' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.009607 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 96.07 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 58.2167 seconds (Warm-up)
Chain 3:                47.0122 seconds (Sampling)
Chain 3:                105.229 seconds (Total)
Chain 3: 
Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 60.1515 seconds (Warm-up)
Chain 2:                47.3321 seconds (Sampling)
Chain 2:                107.484 seconds (Total)
Chain 2: 
Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 69.9765 seconds (Warm-up)
Chain 1:                51.3671 seconds (Sampling)
Chain 1:                121.344 seconds (Total)
Chain 1: 
# stan_hist(fit_post_VFE)
stan_trace(fit_post_VFE)

zP_opt=op_FITC$par[grep("zP",names(op_FITC$par))]
zI_opt=op_FITC$par[grep("zI",names(op_FITC$par))]
data_PI_Z_FITC = data_PI
data_PI_Z_FITC$zP = zP_opt
data_PI_Z_FITC$zI = zI_opt
# plot(sort(zP_opt))
# plot(sort(zI_opt))
zP_opt=op_FITC$par[grep("zP",names(op_FITC$par))]
zI_opt=op_FITC$par[grep("zI",names(op_FITC$par))]
data_PI_Z_FITC = data_PI
data_PI_Z_FITC$zP = zP_opt
data_PI_Z_FITC$zI = zI_opt
# plot(sort(zP_opt))
# plot(sort(zI_opt))
fit_post_FITC=stan(file='STAN/Approximations/FITC/WK2_delta_FITC_fixed_Z.stan',
                  data=data_PI_Z_FITC,
                  chains=3,
                  iter=1000,
                  seed=0
)
starting worker pid=67510 on localhost:11354 at 11:32:59.213
starting worker pid=67524 on localhost:11354 at 11:32:59.413
starting worker pid=67538 on localhost:11354 at 11:32:59.614

SAMPLING FOR MODEL 'WK2_delta_FITC_fixed_Z' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.007932 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 79.32 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'WK2_delta_FITC_fixed_Z' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.007228 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 72.28 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'WK2_delta_FITC_fixed_Z' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.006141 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 61.41 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 60.4329 seconds (Warm-up)
Chain 1:                30.205 seconds (Sampling)
Chain 1:                90.6379 seconds (Total)
Chain 1: 
Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 56.8 seconds (Warm-up)
Chain 3:                38.3528 seconds (Sampling)
Chain 3:                95.1527 seconds (Total)
Chain 3: 
Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 62.991 seconds (Warm-up)
Chain 2:                38.5151 seconds (Sampling)
Chain 2:                101.506 seconds (Total)
Chain 2: 
# stan_hist(fit_post_FITC)
stan_trace(fit_post_FITC)

post_VFE = data.frame(rstan::extract(fit_post_VFE))
post_FITC = data.frame(rstan::extract(fit_post_FITC))
pr = c("R", "C", "sigmaP", "sigmaI")
pVFE = as.vector(as.matrix(post_VFE[,pr]))
pFITC = as.vector(as.matrix(post_FITC[,pr]))
Ns = nrow(post_VFE)
df_plot_post = data.frame(post = c(pVFE, pFITC), par = rep(rep(pr, each = Ns),2), model = c(rep("VFE",length(pr)*Ns), rep("FITC",length(pr)*Ns)))
mod_dat = df_plot_post%>%
  mutate(par = recode(par,
    "R" = "R",
    "C" = "C",
    "sigmaP" = "sigma[P]",
    "sigmaI" = "sigma[Q]"
  ))
df_true_par = data.frame(post=c(Rtrue+Ztrue, Ctrue, 4, 10),par = c("R", "C", "sigma[P]", "sigma[Q]"))
set_lim = data.frame(x=c(0.5,3.1),y=c(600, 600))
df_point_est = data.frame(val=c(op_VFE$par[pr],op_FITC$par[pr]),par = rep(df_true_par$par,2), model = rep(rep(c("VFE","FITC"),each=4),2))
pl_post = ggplot()+
  geom_histogram(data = mod_dat, aes(x=post, fill = par), color="black",bins = 60)+
  facet_grid(model~par, scales = "free", labeller = labeller(par = label_parsed))+
  geom_point(data = set_lim, aes(x=x,y=y), color = "white",alpha=0.000001, size=0.000001)+ # set limits on R and C
  geom_vline(aes(xintercept = post, linetype = "true"), data=df_true_par, size=0.8)+
  geom_vline(aes(xintercept = val, linetype = "map"), data=df_point_est, size=0.8)+
  theme_bw()+ 
  theme(#legend.position = "none",
        legend.title = element_blank(),
        legend.position="bottom",
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        strip.text.x = element_text(size = 13),
        strip.text.y = element_text(size = 13))+ 
  xlab("") + ylab("")+
  scale_fill_manual(
  breaks=c("R", "C", "sigma[P]", "sigma[Q]"),
  values=c("#F8766D","#00BE67","white", "white"),guide = "none")
(pl_post=pl_post + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()))

ggsave("figures/Appr_delta_post.pdf", plot = pl_post, width = 20, height = 12, units = "cm")
nP_pred = 40
ind_P_pred = round(seq(1,length(time),length.out = nP_pred))
tP_pred = time[ind_P_pred]
data_PI_Z_FITC$nP_pred = nP_pred
data_PI_Z_FITC$tP_pred = tP_pred
data_PI_Z_FITC$nI_pred = nP_pred
data_PI_Z_FITC$tI_pred = tP_pred
N_samples = nrow(post_FITC)
data_post_FITC = list(alpha=post_FITC$alpha, rho=post_FITC$rho, alpha_d=post_FITC$alpha_d
                   , rho_d=post_FITC$rho_d, sigmaP=post_FITC$sigmaP, sigmaI=post_FITC$sigmaI
                   , R=post_FITC$R, C=post_FITC$C, N_samples=N_samples
  )
data_pred_FITC = c(data_PI_Z_FITC, data_post_FITC)
pred_FITC = stan(file = 'STAN/Approximations/FITC/FITC_delta_predictions.stan',
                 data = data_pred_FITC,
                 chains = 1, iter = 1, seed=123,
                 algorithm = "Fixed_param")

SAMPLING FOR MODEL 'FITC_delta_predictions' NOW (CHAIN 1).
Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0 seconds (Warm-up)
Chain 1:                3.00098 seconds (Sampling)
Chain 1:                3.00098 seconds (Total)
Chain 1: 
post_mu_CIs.fn = function(post_pred, ci = c(0.05,0.95), time = tP_pred){
  pp = rstan::extract(post_pred)
  dfP = pp$y_P[1,,]
  Psmr = data.frame(mean = colMeans(dfP), 
                    lower = apply(dfP, 2, quantile, probs = ci[1]), 
                    upper = apply(dfP, 2, quantile, probs = ci[2]),
                    time=time)
  dfI = pp$y_I[1,,]
  Ismr = data.frame(mean = colMeans(dfI), 
                    lower = apply(dfI, 2, quantile, probs = ci[1]), 
                    upper = apply(dfI, 2, quantile, probs = ci[2]),
                    time=time)
  return(list(Psmr=Psmr, Ismr=Ismr))
}
data_post_VFE = list(alpha=post_VFE$alpha, rho=post_VFE$rho, alpha_d=post_VFE$alpha_d
                      , rho_d=post_VFE$rho_d, sigmaP=post_VFE$sigmaP, sigmaI=post_VFE$sigmaI
                      , R=post_VFE$R, C=post_VFE$C, N_samples=N_samples
)
data_pred_VFE = c(data_PI_Z_VFE, data_post_VFE)
data_pred_VFE$nP_pred = nP_pred
data_pred_VFE$tP_pred = tP_pred
data_pred_VFE$nI_pred = nP_pred
data_pred_VFE$tI_pred = tP_pred
pred_VFE = stan(file = 'STAN/Approximations/VFE/WK2_delta_VFE_predictions.stan',
                 data = data_pred_VFE,
                 chains = 1, iter = 1, seed=123,
                 algorithm = "Fixed_param")

SAMPLING FOR MODEL 'WK2_delta_VFE_predictions' NOW (CHAIN 1).
Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0 seconds (Warm-up)
Chain 1:                2.47053 seconds (Sampling)
Chain 1:                2.47053 seconds (Total)
Chain 1: 
pp_VFE = rbind(post_mu_CIs.fn(post_pred=pred_VFE)$Psmr, post_mu_CIs.fn(post_pred=pred_VFE)$Ismr)
pp_VFE$output = c(rep("pressure (mmHg)", nP_pred),rep("inflow (ml/min)", nP_pred))
pp_VFE$model = "VFE"
pp_FITC = rbind(post_mu_CIs.fn(post_pred=pred_FITC)$Psmr,post_mu_CIs.fn(post_pred=pred_FITC)$Ismr)
pp_FITC$output = c(rep("pressure (mmHg)", nP_pred),rep("inflow (ml/min)", nP_pred))
pp_FITC$model = "FITC"
pred_df = rbind(pp_VFE, pp_FITC)
# head(pred_df)

Plug in noise estimates

We observe that the VFE model can severely overestimate the noise and therefore result in underfitting. A remedy to this problem is to fix the noise parameter of the functions \(P(t)\) and \(Q(t),\) (\(\sigma_P\) and \(\sigma_I\)). A possible solution for obtaining estimates for the noise parameters is to fit a standard GP model for each dataset \((y_P,t_P)\) and \(y_Q, t_Q\) indepentently and obtain MLE estimates via maximizing the marginal log-likelihood.

df_zP_VFE = data.frame(z=data_PI_Z_VFE$zP, y = rep(60, data_PI_Z_VFE$mP), model = "VFE", output = "pressure (mmHg)")
df_zI_VFE = data.frame(z=data_PI_Z_VFE$zI, y = rep(-100, data_PI_Z_VFE$mI), model = "VFE", output = "inflow (ml/min)")
df_zP_FITC = data.frame(z=data_PI_Z_FITC$zP, y = rep(60, data_PI_Z_FITC$mP), model = "FITC", output = "pressure (mmHg)")
df_zI_FITC = data.frame(z=data_PI_Z_FITC$zI, y = rep(-100, data_PI_Z_FITC$mI), model = "FITC", output = "inflow (ml/min)")
df_z = rbind(df_zP_VFE, df_zI_VFE, df_zP_FITC, df_zI_FITC)
P_true = data.frame(val=Psim, time=time)
P_true$output = "pressure (mmHg)"
I_true = data.frame(val=flow, time=time)
I_true$output = "inflow (ml/min)"
true_out = rbind(P_true, I_true)
obsP = data.frame(value=data_PI$yP, time = data_PI$tP, output = "pressure (mmHg)")
obsI = data.frame(value=data_PI$yI, time = data_PI$tI, output = "inflow (ml/min)")
obs = rbind(obsP,obsI)
pl_pred=ggplot()+
  geom_point(data = obs, aes(y=value, x=time, colour = "observed"), shape = 20)+
  geom_line(data = pred_df, aes(y=mean, x=time, linetype = "mean"), size=0.9)+
  geom_line(data = true_out, aes(y=val, x=time, linetype="true"), size=0.9)+
  geom_ribbon(data = pred_df,aes(ymin=lower, ymax=upper, x=time, fill = "90% CI"), alpha = 0.3)+
  facet_grid(output~model,scales = "free")+
  geom_point(data = df_z, aes(x=z, y=y, shape="Z"), size=3)+
  scale_fill_manual("",values=c("90% CI" = "grey12"))+
  theme_bw()+xlab("time (sec)")+ylab("")+
  scale_shape_manual("", values = c("Z" = 4))+
  theme(#legend.position = "none",
        legend.title = element_blank(),
        legend.position="bottom",
        strip.text.x = element_text(size = 13),
        strip.text.y = element_text(size = 10))
(pl_pred=pl_pred + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()))

ggsave("figures/Appr_delta_pred.pdf", plot = pl_pred, width = 16, height = 10, units = "cm")
nsP = 25
indP = seq(1,data_PI$nP,length.out = nsP)
data_sample_P = list(N=nsP, x = data_PI$tP[indP], y = data_PI$yP[indP])
data_sample_I = list(N=nsP, x = data_PI$tI[indP], y = data_PI$yI[indP])
GP = stan_model('STAN/Approximations/GP_full/GP.stan')
GP_MLE_P=optimizing(GP, data=data_sample_P, hessian=FALSE, verbose=TRUE,seed=0)
Chain 1: Initial log joint probability = -84.3916
Chain 1:     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
Chain 1:       18       -65.767    0.00556334   7.66548e-05           1           1       30   
Chain 1: Optimization terminated normally: 
Chain 1:   Convergence detected: relative gradient magnitude is below tolerance
GP_MLE_P
$par
       rho      alpha      sigma 
 0.2147435 76.0190356  3.3560127 

$value
[1] -65.767

$return_code
[1] 0

$theta_tilde
           rho    alpha    sigma
[1,] 0.2147435 76.01904 3.356013
sigma_P_MLE = GP_MLE_P$par["sigma"] 
GP_MLE_I=optimizing(GP, data=data_sample_I, hessian=FALSE, verbose=TRUE,seed=0)
Chain 1: Initial log joint probability = -656.901
Chain 1:     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
Chain 1:       19      -111.136      0.558199      0.622198           1           1       26   
Chain 1:     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
Chain 1:       39      -111.126      0.321616   0.000908148      0.4931     0.04931       55   
Chain 1:     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
Chain 1:       41      -111.126      0.373032    0.00159772      0.1782           1       58   
Chain 1: Optimization terminated normally: 
Chain 1:   Convergence detected: relative gradient magnitude is below tolerance
GP_MLE_I
$par
         rho        alpha        sigma 
  0.07883712  99.99999925  11.38052032 

$value
[1] -111.1259

$return_code
[1] 0

$theta_tilde
            rho alpha    sigma
[1,] 0.07883712   100 11.38052
sigma_I_MLE = GP_MLE_I$par["sigma"] 
data_pred_VFE_MLE = data_pred_VFE
data_pred_VFE_MLE$sigmaP = sigma_P_MLE
data_pred_VFE_MLE$sigmaI = sigma_I_MLE
pred_VFE_MLE = stan(file = 'STAN/Approximations/VFE/MLE_sigma/WK2_delta_VFE_predictions.stan',
                 data = data_pred_VFE_MLE ,
                 chains = 1, iter = 1, seed=123,
                 algorithm = "Fixed_param")

SAMPLING FOR MODEL 'WK2_delta_VFE_predictions' NOW (CHAIN 1).
Chain 1: Iteration: 1 / 1 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0 seconds (Warm-up)
Chain 1:                2.4886 seconds (Sampling)
Chain 1:                2.4886 seconds (Total)
Chain 1: 
pp_VFE_MLE = rbind(post_mu_CIs.fn(post_pred=pred_VFE_MLE)$Psmr, post_mu_CIs.fn(post_pred=pred_VFE_MLE)$Ismr)
pp_VFE_MLE$output = c(rep("pressure (mmHg)", nP_pred),rep("inflow (ml/min)", nP_pred))
pp_VFE_MLE$model = "VFE fixed noise"
pp_VFE = rbind(post_mu_CIs.fn(post_pred=pred_VFE)$Psmr, post_mu_CIs.fn(post_pred=pred_VFE)$Ismr)
pp_VFE$output = c(rep("pressure (mmHg)", nP_pred),rep("inflow (ml/min)", nP_pred))
pp_VFE$model = "VFE"
pred_df = rbind(pp_VFE_MLE,pp_VFE)
df_zP_VFE = data.frame(z=data_PI_Z_VFE$zP, y = rep(60, data_PI_Z_VFE$mP), model = "VFE", output = "pressure (mmHg)")
df_zI_VFE = data.frame(z=data_PI_Z_VFE$zI, y = rep(-100, data_PI_Z_VFE$mI), model = "VFE", output = "inflow (ml/min)")
df_zP_VFE_MLE = data.frame(z=data_pred_VFE_MLE$zP, y = rep(60, data_pred_VFE_MLE$mP), model = "VFE fixed noise", output = "pressure (mmHg)")
df_zI_VFE_MLE = data.frame(z=data_pred_VFE_MLE$zI, y = rep(-100, data_pred_VFE_MLE$mI), model = "VFE fixed noise", output = "inflow (ml/min)")
df_z = rbind(df_zP_VFE, df_zI_VFE,df_zP_VFE_MLE,df_zI_VFE_MLE)
P_true = data.frame(val=Psim, time=time)
P_true$output = "pressure (mmHg)"
I_true = data.frame(val=flow, time=time)
I_true$output = "inflow (ml/min)"
true_out = rbind(P_true, I_true)
obsP = data.frame(value=data_PI$yP, time = data_PI$tP, output = "pressure (mmHg)")
obsI = data.frame(value=data_PI$yI, time = data_PI$tI, output = "inflow (ml/min)")
obs = rbind(obsP,obsI)
pl_pred=ggplot()+
  geom_point(data = obs, aes(y=value, x=time, colour = "observed"), shape = 20)+
  geom_line(data = pred_df, aes(y=mean, x=time, linetype = "mean"), size=0.9)+
  geom_line(data = true_out, aes(y=val, x=time, linetype="true"), size=0.9)+
  geom_ribbon(data = pred_df,aes(ymin=lower, ymax=upper, x=time, fill = "90% CI"), alpha = 0.3)+
  facet_grid(output~model,scales = "free")+
  geom_point(data = df_z, aes(x=z, y=y, shape="Z"), size=3)+
  scale_fill_manual("",values=c("90% CI" = "grey12"))+
  theme_bw()+xlab("time (sec)")+ylab("")+
  scale_shape_manual("", values = c("Z" = 4))+
  theme(#legend.position = "none",
        legend.title = element_blank(),
        legend.position="bottom",
        strip.text.x = element_text(size = 13),
        strip.text.y = element_text(size = 10))
(pl_pred=pl_pred + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()))

ggsave("figures/Appr_delta_pred_noise.pdf", plot = pl_pred, width = 16, height = 10, units = "cm")
LS0tCnRpdGxlOiAiUzcuMykgQmlnIGRhdGEgYXBwcm94aW1hdGlvbnMgd2l0aCBkaXNjcmVwYW5jeSIKb3V0cHV0OgogIHBkZl9kb2N1bWVudDogZGVmYXVsdAogIGh0bWxfZG9jdW1lbnQ6CiAgICBkZl9wcmludDogcGFnZWQKICBodG1sX25vdGVib29rOiBkZWZhdWx0Ci0tLQoKVGhpcyBub3RlYm9vayBjb250YWlucyB0aGUgY29kZSBvZiB0aGUgcGFwZXIgIkJheWVzaWFuIENhbGlicmF0aW9uIG9mIEltcGVyZmVjdCBDb21wdXRlciBNb2RlbHMgdXNpbmcgUGh5c2ljcy1JbmZvcm1lZCBQcmlvcnMiLiBUaGUgbW9kZWxzIGFyZSBmaXR0ZWQgaW4gcnN0YW4gYW5kIHRoZSBjb2RlIGlzIGF2YWlsYWJsZSBpbiB0aGUgZm9sZGVyICJTVEFOL0FwcHJveGltYXRpb25zIi4gCgojIyMjIExvYWQgcGFja2FnZXMKCmBgYHtyLCBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KAogIG1lc3NhZ2U9RkFMU0UsCiAgd2FybmluZyA9IEZBTFNFLAogIGNvbW1lbnQgPSAnJywgCiAgZmlnLndpZHRoID0gNiwgCiAgZmlnLmhlaWdodCA9IDQsCiAgZmlnLmFsaWduID0gJ2NlbnRlcicKKQpgYGAKCmBgYHtyfQojIHVuY29tbWVudCB0byBpbnN0YWxsCiMgaW5zdGFsbC5wYWNrYWdlcygicnN0YW4iKQojIGluc3RhbGwucGFja2FnZXMoImdncGxvdDIiKQojIGluc3RhbGwucGFja2FnZXMoInRpZHl2ZXJzZSIpCmxpYnJhcnkocnN0YW4pCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeSh0aWR5dmVyc2UpCgpyc3Rhbl9vcHRpb25zKGF1dG9fd3JpdGUgPSBUUlVFKQpvcHRpb25zKG1jLmNvcmVzID0gMykgIyBhbGxvY2F0ZSAzIGNvcmVzIChmb3IgZWFjaCBtb2RlbCB3ZSBydW4gMyBjaGFpbnMgaW4gcGFyYWxsZWwpCiMgTnVtZXJpY2FsIHNpbXVsYXRvciBvZiB0aGUgV0szIG1vZGVsCnNvdXJjZSgiZnVuY3Rpb25zL1dLMmFuZDNfc2ltX2ZuLlIiKQojIExvYWQgZmxvdyBkYXRhIApkID0gcmVhZFJEUygiRGF0YS9JbmZsb3dfdGltZS5yZHMiKQpgYGAKCiMjIyMgUmVhbGl0eSBhbmQgbW9kZWxsaW5nIGNob2ljZSAKClxiZWdpbnthbGlnbn0KICBcbWF0aGNhbHtSfTogXHF1YWQgJiBcZnJhY3tkIFAodCl9e2QgdH0gKyBcZnJhY3tQKHQpfXtSXzJDfSA9IFxmcmFje1EodCl9e0N9IFxsZWZ0ICgxICsgXGZyYWN7IFJfMX17Ul8yfSBccmlnaHQgKSArIFJfMSBcZnJhY3tkIFEodCl9e2R0fSBccXVhZCBcdGV4dHsgKHRoZSBtaXNzcGVzaWZpZWQgbW9kZWwgd2UgdXNlIHRvIGZpdCB0aGUgZGF0YSkgfVx0ZXh0eyBbV0szXSB9XFwKICBcZXRhOiBccXVhZCAgJiAgUSh0KSA9IFxmcmFjezF9e1J9UCh0KSArIEMgXGZyYWN7ZFAodCl9e2R0fSBccXVhZCAgXHRleHR7ICh0aGUgbW9kZWwgd2UgdXNlIHRvIHNpbXVsYXRlIGRhdGEpIH0gXHRleHR7IFtXSzJdIH0KXGVuZHthbGlnbn0KCgoKYGBge3IgZXZhbD1UUlVFfQpSdHJ1ZSA9IDE7IEN0cnVlID0gMS4xOyBadHJ1ZSA9IDAuMDUgCmZsb3cgPSBkJGluZmxvdyowLjk1CnRpbWUgPSBkJHRpbWUKblAgPSA5MCAjIG51bWJlciBvZiBwcmVzc3VyZSBkYXRhCm5JID0gMTAwIyBudW1iZXIgb2YgaW5mbG93IGRhdGEKbmMgPSAxICAjIG51bWJlciBvZiBjYXJkaWFjIGN5Y2xlcwpuZmxvdyA9IGxlbmd0aChmbG93KQojIDEuIHNpbXVsYXRlIFdLMyBkYXRhIChSPVJfMiwgWj1SXzEpClBzaW0gPSBXSzNfc2ltdWxhdGUoZmxvdyA9IGZsb3csIHRpbWUgPSB0aW1lLCBSID0gUnRydWUsIEMgPSBDdHJ1ZSwgWj1adHJ1ZSkgIyBzaW11bGF0ZSBXSzMgZGF0YSBmb3IgYSBnaXZlbiBmbG93IFEodCkKUF90cnVlID0gUHNpbQojIDIuIGNob29zZSBwcmVzc3VyZSBhbmQgaW5mbG93IGluZGljZXMKaW5kUCA9IHJvdW5kKHNlcSgxLCBuZmxvdywgbGVuZ3RoLm91dCA9IG5QKSk7IGluZEkgPSByb3VuZChzZXEoMSwgbmZsb3csIGxlbmd0aC5vdXQgPSBuSSkpCnlQX3JlYWwgPSBQc2ltW2luZFBdOyB5SV9yZWFsID0gZmxvd1tpbmRJXSAjIG5vaXNlIGZyZWUgZmltdWxhdGVkIHByZXNzdXJlIGFuZCBmbG93CiMgMy4gQWRkIG5vaXNlCiMgc2V0LnNlZWQoMCkKc2V0LnNlZWQoMSkKUG5vaXNlID0gcm5vcm0oblAqbmMsIDAsIDQpICMgc2FtcGxlIHByZXNzdXJlIG5vaXNlIGZyb20gTigwLCA0XjIpCklub2lzZSA9IHJub3JtKG5JKm5jLCAwLCAxMCkgIyBzYW1wbGUgZmxvdyBub2lzZSBmcm9tIE4oMCwxMF4yKQp5UF9yZWFsID0gcmVwKHlQX3JlYWwsbmMpIAp5SV9yZWFsID0gcmVwKHlJX3JlYWwsbmMpCiMgNC4gc3RvcmUgaW5kaXZpZHVhbCBkYXRhIGluIHRoZSBwb3B1bGF0aW9uIG1hdHJpY2VzCnlQID0geVBfcmVhbCArIFBub2lzZSAjIGFkZCBub2lzZQp5SSA9IHlJX3JlYWwgKyBJbm9pc2UgIyBhZGQgbm9pc2UKdFAgPSB0aW1lW2luZFBdICMgY29ycmVzcG9uZGluZyB0aW1lIChzeW5jaHJvbml6ZWQgZm9yIHRoZSB0d28gY3ljbGVzKQp0SSA9IHRpbWVbaW5kSV0gIyBjb3JyZXNwb25kaW5nIHRpbWUgKHN5bmNocm9uaXplZCBmb3IgdGhlIHR3byBjeWNsZXMpCgpkYXRhX1BJID0gbGlzdChuUD1uYypuUCwgbkk9bmMqbkksIHRQPXJlcCh0UCxuYyksIHRJPXJlcCh0SSxuYyksIHlQPXlQLCB5ST15SSwgbVA9MTIsIG1JPTEwKQpgYGAKCgpgYGB7cn0KV0syX1ZGRSA9IHN0YW5fbW9kZWwoJ1NUQU4vQXBwcm94aW1hdGlvbnMvVkZFL1dLMl9kZWx0YV9WRkUuc3RhbicpCmtwID0ga21lYW5zKGRhdGEuZnJhbWUoeD1kYXRhX1BJJHRQKSwgY2VudGVycyA9IGRhdGFfUEkkbVApCmtpID0ga21lYW5zKGRhdGEuZnJhbWUoeD1kYXRhX1BJJHRJKSwgY2VudGVycyA9IGRhdGFfUEkkbUkpCgppbml0ID0gbGlzdCgielAiID0gYXMudmVjdG9yKGtwJGNlbnRlcnMpLCAiekkiID0gYXMudmVjdG9yKGtpJGNlbnRlcnMpKQpvcF9WRkU9b3B0aW1pemluZyhXSzJfVkZFLCBkYXRhPWRhdGFfUEksIGhlc3NpYW49RkFMU0UsIGluaXQgPSBpbml0LCB2ZXJib3NlPVRSVUUsIHNlZWQ9MCkKb3BfVkZFCmBgYAoKYGBge3J9CnpQX29wdF9WRkU9b3BfVkZFJHBhcltncmVwKCJ6UCIsbmFtZXMob3BfVkZFJHBhcikpXQp6SV9vcHRfVkZFPW9wX1ZGRSRwYXJbZ3JlcCgiekkiLG5hbWVzKG9wX1ZGRSRwYXIpKV0KIyBwbG90KHNvcnQoelBfb3B0X1ZGRSkpCiMgcGxvdChzb3J0KHpJX29wdF9WRkUpKQpkYXRhX1BJX1pfVkZFPSBkYXRhX1BJCmRhdGFfUElfWl9WRkUkelAgPSB6UF9vcHRfVkZFCmRhdGFfUElfWl9WRkUkekkgPSB6SV9vcHRfVkZFCgpmaXRfcG9zdF9WRkU9c3RhbihmaWxlPSdTVEFOL0FwcHJveGltYXRpb25zL1ZGRS9XSzJfZGVsdGFfVkZFX2ZpeGVkX1ouc3RhbicsCiAgICAgICAgICAgICAgICAgIGRhdGE9ZGF0YV9QSV9aX1ZGRSwKICAgICAgICAgICAgICAgICAgY2hhaW5zPTMsCiAgICAgICAgICAgICAgICAgIGl0ZXI9MTAwMCwKICAgICAgICAgICAgICAgICAgc2VlZD0wCikKIyBzdGFuX2hpc3QoZml0X3Bvc3RfVkZFKQpzdGFuX3RyYWNlKGZpdF9wb3N0X1ZGRSkKYGBgCgojIyMgRklUQyAKCmBgYHtyfQpXSzJfRklUQyA9IHN0YW5fbW9kZWwoJ1NUQU4vQXBwcm94aW1hdGlvbnMvRklUQy9XSzJfZGVsdGFfRklUQy5zdGFuJykKb3BfRklUQz1vcHRpbWl6aW5nKFdLMl9GSVRDLCBkYXRhPWRhdGFfUEksIGhlc3NpYW49RkFMU0UsIHZlcmJvc2U9VFJVRSxpbml0PWluaXQsc2VlZD0zMSkKb3BfRklUQwpgYGAKCmBgYHtyfQp6UF9vcHQ9b3BfRklUQyRwYXJbZ3JlcCgielAiLG5hbWVzKG9wX0ZJVEMkcGFyKSldCnpJX29wdD1vcF9GSVRDJHBhcltncmVwKCJ6SSIsbmFtZXMob3BfRklUQyRwYXIpKV0KZGF0YV9QSV9aX0ZJVEMgPSBkYXRhX1BJCmRhdGFfUElfWl9GSVRDJHpQID0gelBfb3B0CmRhdGFfUElfWl9GSVRDJHpJID0geklfb3B0CiMgcGxvdChzb3J0KHpQX29wdCkpCiMgcGxvdChzb3J0KHpJX29wdCkpCmBgYAoKCmBgYHtyfQpmaXRfcG9zdF9GSVRDPXN0YW4oZmlsZT0nU1RBTi9BcHByb3hpbWF0aW9ucy9GSVRDL1dLMl9kZWx0YV9GSVRDX2ZpeGVkX1ouc3RhbicsCiAgICAgICAgICAgICAgICAgIGRhdGE9ZGF0YV9QSV9aX0ZJVEMsCiAgICAgICAgICAgICAgICAgIGNoYWlucz0zLAogICAgICAgICAgICAgICAgICBpdGVyPTEwMDAsCiAgICAgICAgICAgICAgICAgIHNlZWQ9MAopCgojIHN0YW5faGlzdChmaXRfcG9zdF9GSVRDKQpzdGFuX3RyYWNlKGZpdF9wb3N0X0ZJVEMpCmBgYAoKCgoKCgoKCgpgYGB7cn0KcG9zdF9WRkUgPSBkYXRhLmZyYW1lKHJzdGFuOjpleHRyYWN0KGZpdF9wb3N0X1ZGRSkpCnBvc3RfRklUQyA9IGRhdGEuZnJhbWUocnN0YW46OmV4dHJhY3QoZml0X3Bvc3RfRklUQykpCgpwciA9IGMoIlIiLCAiQyIsICJzaWdtYVAiLCAic2lnbWFJIikKcFZGRSA9IGFzLnZlY3Rvcihhcy5tYXRyaXgocG9zdF9WRkVbLHByXSkpCnBGSVRDID0gYXMudmVjdG9yKGFzLm1hdHJpeChwb3N0X0ZJVENbLHByXSkpCk5zID0gbnJvdyhwb3N0X1ZGRSkKCmRmX3Bsb3RfcG9zdCA9IGRhdGEuZnJhbWUocG9zdCA9IGMocFZGRSwgcEZJVEMpLCBwYXIgPSByZXAocmVwKHByLCBlYWNoID0gTnMpLDIpLCBtb2RlbCA9IGMocmVwKCJWRkUiLGxlbmd0aChwcikqTnMpLCByZXAoIkZJVEMiLGxlbmd0aChwcikqTnMpKSkKYGBgCgpgYGB7cixmaWcud2lkdGg9MTAsIGZpZy5oZWlnaHQ9NX0KbW9kX2RhdCA9IGRmX3Bsb3RfcG9zdCU+JQogIG11dGF0ZShwYXIgPSByZWNvZGUocGFyLAogICAgIlIiID0gIlIiLAogICAgIkMiID0gIkMiLAogICAgInNpZ21hUCIgPSAic2lnbWFbUF0iLAogICAgInNpZ21hSSIgPSAic2lnbWFbUV0iCiAgKSkKZGZfdHJ1ZV9wYXIgPSBkYXRhLmZyYW1lKHBvc3Q9YyhSdHJ1ZStadHJ1ZSwgQ3RydWUsIDQsIDEwKSxwYXIgPSBjKCJSIiwgIkMiLCAic2lnbWFbUF0iLCAic2lnbWFbUV0iKSkKc2V0X2xpbSA9IGRhdGEuZnJhbWUoeD1jKDAuNSwzLjEpLHk9Yyg2MDAsIDYwMCkpCmRmX3BvaW50X2VzdCA9IGRhdGEuZnJhbWUodmFsPWMob3BfVkZFJHBhcltwcl0sb3BfRklUQyRwYXJbcHJdKSxwYXIgPSByZXAoZGZfdHJ1ZV9wYXIkcGFyLDIpLCBtb2RlbCA9IHJlcChyZXAoYygiVkZFIiwiRklUQyIpLGVhY2g9NCksMikpCnBsX3Bvc3QgPSBnZ3Bsb3QoKSsKICBnZW9tX2hpc3RvZ3JhbShkYXRhID0gbW9kX2RhdCwgYWVzKHg9cG9zdCwgZmlsbCA9IHBhciksIGNvbG9yPSJibGFjayIsYmlucyA9IDYwKSsKICBmYWNldF9ncmlkKG1vZGVsfnBhciwgc2NhbGVzID0gImZyZWUiLCBsYWJlbGxlciA9IGxhYmVsbGVyKHBhciA9IGxhYmVsX3BhcnNlZCkpKwogIGdlb21fcG9pbnQoZGF0YSA9IHNldF9saW0sIGFlcyh4PXgseT15KSwgY29sb3IgPSAid2hpdGUiLGFscGhhPTAuMDAwMDAxLCBzaXplPTAuMDAwMDAxKSsgIyBzZXQgbGltaXRzIG9uIFIgYW5kIEMKICBnZW9tX3ZsaW5lKGFlcyh4aW50ZXJjZXB0ID0gcG9zdCwgbGluZXR5cGUgPSAidHJ1ZSIpLCBkYXRhPWRmX3RydWVfcGFyLCBzaXplPTAuOCkrCiAgZ2VvbV92bGluZShhZXMoeGludGVyY2VwdCA9IHZhbCwgbGluZXR5cGUgPSAibWFwIiksIGRhdGE9ZGZfcG9pbnRfZXN0LCBzaXplPTAuOCkrCiAgdGhlbWVfYncoKSsgCiAgdGhlbWUoI2xlZ2VuZC5wb3NpdGlvbiA9ICJub25lIiwKICAgICAgICBsZWdlbmQudGl0bGUgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgbGVnZW5kLnBvc2l0aW9uPSJib3R0b20iLAogICAgICAgIGF4aXMudGl0bGUueT1lbGVtZW50X2JsYW5rKCksCiAgICAgICAgYXhpcy50ZXh0Lnk9ZWxlbWVudF9ibGFuaygpLAogICAgICAgIGF4aXMudGlja3MueT1lbGVtZW50X2JsYW5rKCksCiAgICAgICAgc3RyaXAudGV4dC54ID0gZWxlbWVudF90ZXh0KHNpemUgPSAxMyksCiAgICAgICAgc3RyaXAudGV4dC55ID0gZWxlbWVudF90ZXh0KHNpemUgPSAxMykpKyAKICB4bGFiKCIiKSArIHlsYWIoIiIpKwogIHNjYWxlX2ZpbGxfbWFudWFsKAogIGJyZWFrcz1jKCJSIiwgIkMiLCAic2lnbWFbUF0iLCAic2lnbWFbUV0iKSwKICB2YWx1ZXM9YygiI0Y4NzY2RCIsIiMwMEJFNjciLCJ3aGl0ZSIsICJ3aGl0ZSIpLGd1aWRlID0gIm5vbmUiKQoocGxfcG9zdD1wbF9wb3N0ICsgdGhlbWUocGFuZWwuZ3JpZC5tYWpvciA9IGVsZW1lbnRfYmxhbmsoKSwgcGFuZWwuZ3JpZC5taW5vciA9IGVsZW1lbnRfYmxhbmsoKSkpCmdnc2F2ZSgiZmlndXJlcy9BcHByX2RlbHRhX3Bvc3QucGRmIiwgcGxvdCA9IHBsX3Bvc3QsIHdpZHRoID0gMjAsIGhlaWdodCA9IDEyLCB1bml0cyA9ICJjbSIpCmBgYAoKYGBge3J9Cm5QX3ByZWQgPSA0MAppbmRfUF9wcmVkID0gcm91bmQoc2VxKDEsbGVuZ3RoKHRpbWUpLGxlbmd0aC5vdXQgPSBuUF9wcmVkKSkKdFBfcHJlZCA9IHRpbWVbaW5kX1BfcHJlZF0KZGF0YV9QSV9aX0ZJVEMkblBfcHJlZCA9IG5QX3ByZWQKZGF0YV9QSV9aX0ZJVEMkdFBfcHJlZCA9IHRQX3ByZWQKZGF0YV9QSV9aX0ZJVEMkbklfcHJlZCA9IG5QX3ByZWQKZGF0YV9QSV9aX0ZJVEMkdElfcHJlZCA9IHRQX3ByZWQKCk5fc2FtcGxlcyA9IG5yb3cocG9zdF9GSVRDKQpkYXRhX3Bvc3RfRklUQyA9IGxpc3QoYWxwaGE9cG9zdF9GSVRDJGFscGhhLCByaG89cG9zdF9GSVRDJHJobywgYWxwaGFfZD1wb3N0X0ZJVEMkYWxwaGFfZAogICAgICAgICAgICAgICAgICAgLCByaG9fZD1wb3N0X0ZJVEMkcmhvX2QsIHNpZ21hUD1wb3N0X0ZJVEMkc2lnbWFQLCBzaWdtYUk9cG9zdF9GSVRDJHNpZ21hSQogICAgICAgICAgICAgICAgICAgLCBSPXBvc3RfRklUQyRSLCBDPXBvc3RfRklUQyRDLCBOX3NhbXBsZXM9Tl9zYW1wbGVzCiAgKQoKZGF0YV9wcmVkX0ZJVEMgPSBjKGRhdGFfUElfWl9GSVRDLCBkYXRhX3Bvc3RfRklUQykKYGBgCgpgYGB7cn0KcHJlZF9GSVRDID0gc3RhbihmaWxlID0gJ1NUQU4vQXBwcm94aW1hdGlvbnMvRklUQy9GSVRDX2RlbHRhX3ByZWRpY3Rpb25zLnN0YW4nLAogICAgICAgICAgICAgICAgIGRhdGEgPSBkYXRhX3ByZWRfRklUQywKICAgICAgICAgICAgICAgICBjaGFpbnMgPSAxLCBpdGVyID0gMSwgc2VlZD0xMjMsCiAgICAgICAgICAgICAgICAgYWxnb3JpdGhtID0gIkZpeGVkX3BhcmFtIikKYGBgCgpgYGB7cn0KcG9zdF9tdV9DSXMuZm4gPSBmdW5jdGlvbihwb3N0X3ByZWQsIGNpID0gYygwLjA1LDAuOTUpLCB0aW1lID0gdFBfcHJlZCl7CiAgcHAgPSByc3Rhbjo6ZXh0cmFjdChwb3N0X3ByZWQpCiAgZGZQID0gcHAkeV9QWzEsLF0KICBQc21yID0gZGF0YS5mcmFtZShtZWFuID0gY29sTWVhbnMoZGZQKSwgCiAgICAgICAgICAgICAgICAgICAgbG93ZXIgPSBhcHBseShkZlAsIDIsIHF1YW50aWxlLCBwcm9icyA9IGNpWzFdKSwgCiAgICAgICAgICAgICAgICAgICAgdXBwZXIgPSBhcHBseShkZlAsIDIsIHF1YW50aWxlLCBwcm9icyA9IGNpWzJdKSwKICAgICAgICAgICAgICAgICAgICB0aW1lPXRpbWUpCiAgZGZJID0gcHAkeV9JWzEsLF0KICBJc21yID0gZGF0YS5mcmFtZShtZWFuID0gY29sTWVhbnMoZGZJKSwgCiAgICAgICAgICAgICAgICAgICAgbG93ZXIgPSBhcHBseShkZkksIDIsIHF1YW50aWxlLCBwcm9icyA9IGNpWzFdKSwgCiAgICAgICAgICAgICAgICAgICAgdXBwZXIgPSBhcHBseShkZkksIDIsIHF1YW50aWxlLCBwcm9icyA9IGNpWzJdKSwKICAgICAgICAgICAgICAgICAgICB0aW1lPXRpbWUpCiAgcmV0dXJuKGxpc3QoUHNtcj1Qc21yLCBJc21yPUlzbXIpKQp9CmBgYAoKYGBge3J9CmRhdGFfcG9zdF9WRkUgPSBsaXN0KGFscGhhPXBvc3RfVkZFJGFscGhhLCByaG89cG9zdF9WRkUkcmhvLCBhbHBoYV9kPXBvc3RfVkZFJGFscGhhX2QKICAgICAgICAgICAgICAgICAgICAgICwgcmhvX2Q9cG9zdF9WRkUkcmhvX2QsIHNpZ21hUD1wb3N0X1ZGRSRzaWdtYVAsIHNpZ21hST1wb3N0X1ZGRSRzaWdtYUkKICAgICAgICAgICAgICAgICAgICAgICwgUj1wb3N0X1ZGRSRSLCBDPXBvc3RfVkZFJEMsIE5fc2FtcGxlcz1OX3NhbXBsZXMKKQoKCmRhdGFfcHJlZF9WRkUgPSBjKGRhdGFfUElfWl9WRkUsIGRhdGFfcG9zdF9WRkUpCmRhdGFfcHJlZF9WRkUkblBfcHJlZCA9IG5QX3ByZWQKZGF0YV9wcmVkX1ZGRSR0UF9wcmVkID0gdFBfcHJlZApkYXRhX3ByZWRfVkZFJG5JX3ByZWQgPSBuUF9wcmVkCmRhdGFfcHJlZF9WRkUkdElfcHJlZCA9IHRQX3ByZWQKCmBgYAoKYGBge3J9CnByZWRfVkZFID0gc3RhbihmaWxlID0gJ1NUQU4vQXBwcm94aW1hdGlvbnMvVkZFL1dLMl9kZWx0YV9WRkVfcHJlZGljdGlvbnMuc3RhbicsCiAgICAgICAgICAgICAgICAgZGF0YSA9IGRhdGFfcHJlZF9WRkUsCiAgICAgICAgICAgICAgICAgY2hhaW5zID0gMSwgaXRlciA9IDEsIHNlZWQ9MTIzLAogICAgICAgICAgICAgICAgIGFsZ29yaXRobSA9ICJGaXhlZF9wYXJhbSIpCmBgYAoKYGBge3J9CnBwX1ZGRSA9IHJiaW5kKHBvc3RfbXVfQ0lzLmZuKHBvc3RfcHJlZD1wcmVkX1ZGRSkkUHNtciwgcG9zdF9tdV9DSXMuZm4ocG9zdF9wcmVkPXByZWRfVkZFKSRJc21yKQpwcF9WRkUkb3V0cHV0ID0gYyhyZXAoInByZXNzdXJlIChtbUhnKSIsIG5QX3ByZWQpLHJlcCgiaW5mbG93IChtbC9taW4pIiwgblBfcHJlZCkpCnBwX1ZGRSRtb2RlbCA9ICJWRkUiCnBwX0ZJVEMgPSByYmluZChwb3N0X211X0NJcy5mbihwb3N0X3ByZWQ9cHJlZF9GSVRDKSRQc21yLHBvc3RfbXVfQ0lzLmZuKHBvc3RfcHJlZD1wcmVkX0ZJVEMpJElzbXIpCnBwX0ZJVEMkb3V0cHV0ID0gYyhyZXAoInByZXNzdXJlIChtbUhnKSIsIG5QX3ByZWQpLHJlcCgiaW5mbG93IChtbC9taW4pIiwgblBfcHJlZCkpCnBwX0ZJVEMkbW9kZWwgPSAiRklUQyIKcHJlZF9kZiA9IHJiaW5kKHBwX1ZGRSwgcHBfRklUQykKIyBoZWFkKHByZWRfZGYpCmBgYAoKCgpgYGB7cixmaWcud2lkdGg9OH0KZGZfelBfVkZFID0gZGF0YS5mcmFtZSh6PWRhdGFfUElfWl9WRkUkelAsIHkgPSByZXAoNjAsIGRhdGFfUElfWl9WRkUkbVApLCBtb2RlbCA9ICJWRkUiLCBvdXRwdXQgPSAicHJlc3N1cmUgKG1tSGcpIikKZGZfeklfVkZFID0gZGF0YS5mcmFtZSh6PWRhdGFfUElfWl9WRkUkekksIHkgPSByZXAoLTEwMCwgZGF0YV9QSV9aX1ZGRSRtSSksIG1vZGVsID0gIlZGRSIsIG91dHB1dCA9ICJpbmZsb3cgKG1sL21pbikiKQpkZl96UF9GSVRDID0gZGF0YS5mcmFtZSh6PWRhdGFfUElfWl9GSVRDJHpQLCB5ID0gcmVwKDYwLCBkYXRhX1BJX1pfRklUQyRtUCksIG1vZGVsID0gIkZJVEMiLCBvdXRwdXQgPSAicHJlc3N1cmUgKG1tSGcpIikKZGZfeklfRklUQyA9IGRhdGEuZnJhbWUoej1kYXRhX1BJX1pfRklUQyR6SSwgeSA9IHJlcCgtMTAwLCBkYXRhX1BJX1pfRklUQyRtSSksIG1vZGVsID0gIkZJVEMiLCBvdXRwdXQgPSAiaW5mbG93IChtbC9taW4pIikKZGZfeiA9IHJiaW5kKGRmX3pQX1ZGRSwgZGZfeklfVkZFLCBkZl96UF9GSVRDLCBkZl96SV9GSVRDKQpQX3RydWUgPSBkYXRhLmZyYW1lKHZhbD1Qc2ltLCB0aW1lPXRpbWUpClBfdHJ1ZSRvdXRwdXQgPSAicHJlc3N1cmUgKG1tSGcpIgpJX3RydWUgPSBkYXRhLmZyYW1lKHZhbD1mbG93LCB0aW1lPXRpbWUpCklfdHJ1ZSRvdXRwdXQgPSAiaW5mbG93IChtbC9taW4pIgp0cnVlX291dCA9IHJiaW5kKFBfdHJ1ZSwgSV90cnVlKQoKb2JzUCA9IGRhdGEuZnJhbWUodmFsdWU9ZGF0YV9QSSR5UCwgdGltZSA9IGRhdGFfUEkkdFAsIG91dHB1dCA9ICJwcmVzc3VyZSAobW1IZykiKQpvYnNJID0gZGF0YS5mcmFtZSh2YWx1ZT1kYXRhX1BJJHlJLCB0aW1lID0gZGF0YV9QSSR0SSwgb3V0cHV0ID0gImluZmxvdyAobWwvbWluKSIpCm9icyA9IHJiaW5kKG9ic1Asb2JzSSkKCnBsX3ByZWQ9Z2dwbG90KCkrCiAgZ2VvbV9wb2ludChkYXRhID0gb2JzLCBhZXMoeT12YWx1ZSwgeD10aW1lLCBjb2xvdXIgPSAib2JzZXJ2ZWQiKSwgc2hhcGUgPSAyMCkrCiAgZ2VvbV9saW5lKGRhdGEgPSBwcmVkX2RmLCBhZXMoeT1tZWFuLCB4PXRpbWUsIGxpbmV0eXBlID0gIm1lYW4iKSwgc2l6ZT0wLjkpKwogIGdlb21fbGluZShkYXRhID0gdHJ1ZV9vdXQsIGFlcyh5PXZhbCwgeD10aW1lLCBsaW5ldHlwZT0idHJ1ZSIpLCBzaXplPTAuOSkrCiAgZ2VvbV9yaWJib24oZGF0YSA9IHByZWRfZGYsYWVzKHltaW49bG93ZXIsIHltYXg9dXBwZXIsIHg9dGltZSwgZmlsbCA9ICI5MCUgQ0kiKSwgYWxwaGEgPSAwLjMpKwogIGZhY2V0X2dyaWQob3V0cHV0fm1vZGVsLHNjYWxlcyA9ICJmcmVlIikrCiAgZ2VvbV9wb2ludChkYXRhID0gZGZfeiwgYWVzKHg9eiwgeT15LCBzaGFwZT0iWiIpLCBzaXplPTMpKwogIHNjYWxlX2ZpbGxfbWFudWFsKCIiLHZhbHVlcz1jKCI5MCUgQ0kiID0gImdyZXkxMiIpKSsKICB0aGVtZV9idygpK3hsYWIoInRpbWUgKHNlYykiKSt5bGFiKCIiKSsKICBzY2FsZV9zaGFwZV9tYW51YWwoIiIsIHZhbHVlcyA9IGMoIloiID0gNCkpKwogIHRoZW1lKCNsZWdlbmQucG9zaXRpb24gPSAibm9uZSIsCiAgICAgICAgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgIGxlZ2VuZC5wb3NpdGlvbj0iYm90dG9tIiwKICAgICAgICBzdHJpcC50ZXh0LnggPSBlbGVtZW50X3RleHQoc2l6ZSA9IDEzKSwKICAgICAgICBzdHJpcC50ZXh0LnkgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDEwKSkKKHBsX3ByZWQ9cGxfcHJlZCArIHRoZW1lKHBhbmVsLmdyaWQubWFqb3IgPSBlbGVtZW50X2JsYW5rKCksIHBhbmVsLmdyaWQubWlub3IgPSBlbGVtZW50X2JsYW5rKCkpKQoKZ2dzYXZlKCJmaWd1cmVzL0FwcHJfZGVsdGFfcHJlZC5wZGYiLCBwbG90ID0gcGxfcHJlZCwgd2lkdGggPSAxNiwgaGVpZ2h0ID0gMTAsIHVuaXRzID0gImNtIikKYGBgCgojIyMgUGx1ZyBpbiBub2lzZSAgZXN0aW1hdGVzIAoKV2Ugb2JzZXJ2ZSB0aGF0IHRoZSBWRkUgbW9kZWwgY2FuIHNldmVyZWx5IG92ZXJlc3RpbWF0ZSB0aGUgbm9pc2UgYW5kIHRoZXJlZm9yZSByZXN1bHQgaW4gdW5kZXJmaXR0aW5nLiBBIHJlbWVkeSB0byB0aGlzIHByb2JsZW0gaXMgdG8gZml4IHRoZSBub2lzZSBwYXJhbWV0ZXIgb2YgdGhlIGZ1bmN0aW9ucyAkUCh0KSQgYW5kICRRKHQpLCQgKCRcc2lnbWFfUCQgYW5kICRcc2lnbWFfSSQpLiBBIHBvc3NpYmxlIHNvbHV0aW9uIGZvciBvYnRhaW5pbmcgZXN0aW1hdGVzIGZvciB0aGUgbm9pc2UgcGFyYW1ldGVycyBpcyB0byBmaXQgYSBzdGFuZGFyZCBHUCBtb2RlbCBmb3IgZWFjaCBkYXRhc2V0ICQoeV9QLHRfUCkkIGFuZCAkeV9RLCB0X1EkIGluZGVwZW50ZW50bHkgYW5kIG9idGFpbiBNTEUgZXN0aW1hdGVzIHZpYSBtYXhpbWl6aW5nIHRoZSBtYXJnaW5hbCBsb2ctbGlrZWxpaG9vZC4KCgpgYGB7cn0KbnNQID0gMjUKaW5kUCA9IHNlcSgxLGRhdGFfUEkkblAsbGVuZ3RoLm91dCA9IG5zUCkKZGF0YV9zYW1wbGVfUCA9IGxpc3QoTj1uc1AsIHggPSBkYXRhX1BJJHRQW2luZFBdLCB5ID0gZGF0YV9QSSR5UFtpbmRQXSkKZGF0YV9zYW1wbGVfSSA9IGxpc3QoTj1uc1AsIHggPSBkYXRhX1BJJHRJW2luZFBdLCB5ID0gZGF0YV9QSSR5SVtpbmRQXSkKCkdQID0gc3Rhbl9tb2RlbCgnU1RBTi9BcHByb3hpbWF0aW9ucy9HUF9mdWxsL0dQLnN0YW4nKQoKR1BfTUxFX1A9b3B0aW1pemluZyhHUCwgZGF0YT1kYXRhX3NhbXBsZV9QLCBoZXNzaWFuPUZBTFNFLCB2ZXJib3NlPVRSVUUsc2VlZD0wKQpHUF9NTEVfUAoKc2lnbWFfUF9NTEUgPSBHUF9NTEVfUCRwYXJbInNpZ21hIl0gCgpHUF9NTEVfST1vcHRpbWl6aW5nKEdQLCBkYXRhPWRhdGFfc2FtcGxlX0ksIGhlc3NpYW49RkFMU0UsIHZlcmJvc2U9VFJVRSxzZWVkPTApCkdQX01MRV9JCgpzaWdtYV9JX01MRSA9IEdQX01MRV9JJHBhclsic2lnbWEiXSAKCmBgYApgYGB7cn0KZGF0YV9wcmVkX1ZGRV9NTEUgPSBkYXRhX3ByZWRfVkZFCmRhdGFfcHJlZF9WRkVfTUxFJHNpZ21hUCA9IHNpZ21hX1BfTUxFCmRhdGFfcHJlZF9WRkVfTUxFJHNpZ21hSSA9IHNpZ21hX0lfTUxFCnByZWRfVkZFX01MRSA9IHN0YW4oZmlsZSA9ICdTVEFOL0FwcHJveGltYXRpb25zL1ZGRS9NTEVfc2lnbWEvV0syX2RlbHRhX1ZGRV9wcmVkaWN0aW9ucy5zdGFuJywKICAgICAgICAgICAgICAgICBkYXRhID0gZGF0YV9wcmVkX1ZGRV9NTEUgLAogICAgICAgICAgICAgICAgIGNoYWlucyA9IDEsIGl0ZXIgPSAxLCBzZWVkPTEyMywKICAgICAgICAgICAgICAgICBhbGdvcml0aG0gPSAiRml4ZWRfcGFyYW0iKQpgYGAKCgpgYGB7cn0KcHBfVkZFX01MRSA9IHJiaW5kKHBvc3RfbXVfQ0lzLmZuKHBvc3RfcHJlZD1wcmVkX1ZGRV9NTEUpJFBzbXIsIHBvc3RfbXVfQ0lzLmZuKHBvc3RfcHJlZD1wcmVkX1ZGRV9NTEUpJElzbXIpCnBwX1ZGRV9NTEUkb3V0cHV0ID0gYyhyZXAoInByZXNzdXJlIChtbUhnKSIsIG5QX3ByZWQpLHJlcCgiaW5mbG93IChtbC9taW4pIiwgblBfcHJlZCkpCnBwX1ZGRV9NTEUkbW9kZWwgPSAiVkZFIGZpeGVkIG5vaXNlIgpwcF9WRkUgPSByYmluZChwb3N0X211X0NJcy5mbihwb3N0X3ByZWQ9cHJlZF9WRkUpJFBzbXIsIHBvc3RfbXVfQ0lzLmZuKHBvc3RfcHJlZD1wcmVkX1ZGRSkkSXNtcikKcHBfVkZFJG91dHB1dCA9IGMocmVwKCJwcmVzc3VyZSAobW1IZykiLCBuUF9wcmVkKSxyZXAoImluZmxvdyAobWwvbWluKSIsIG5QX3ByZWQpKQpwcF9WRkUkbW9kZWwgPSAiVkZFIgpwcmVkX2RmID0gcmJpbmQocHBfVkZFX01MRSxwcF9WRkUpCmRmX3pQX1ZGRSA9IGRhdGEuZnJhbWUoej1kYXRhX1BJX1pfVkZFJHpQLCB5ID0gcmVwKDYwLCBkYXRhX1BJX1pfVkZFJG1QKSwgbW9kZWwgPSAiVkZFIiwgb3V0cHV0ID0gInByZXNzdXJlIChtbUhnKSIpCmRmX3pJX1ZGRSA9IGRhdGEuZnJhbWUoej1kYXRhX1BJX1pfVkZFJHpJLCB5ID0gcmVwKC0xMDAsIGRhdGFfUElfWl9WRkUkbUkpLCBtb2RlbCA9ICJWRkUiLCBvdXRwdXQgPSAiaW5mbG93IChtbC9taW4pIikKZGZfelBfVkZFX01MRSA9IGRhdGEuZnJhbWUoej1kYXRhX3ByZWRfVkZFX01MRSR6UCwgeSA9IHJlcCg2MCwgZGF0YV9wcmVkX1ZGRV9NTEUkbVApLCBtb2RlbCA9ICJWRkUgZml4ZWQgbm9pc2UiLCBvdXRwdXQgPSAicHJlc3N1cmUgKG1tSGcpIikKZGZfeklfVkZFX01MRSA9IGRhdGEuZnJhbWUoej1kYXRhX3ByZWRfVkZFX01MRSR6SSwgeSA9IHJlcCgtMTAwLCBkYXRhX3ByZWRfVkZFX01MRSRtSSksIG1vZGVsID0gIlZGRSBmaXhlZCBub2lzZSIsIG91dHB1dCA9ICJpbmZsb3cgKG1sL21pbikiKQpkZl96ID0gcmJpbmQoZGZfelBfVkZFLCBkZl96SV9WRkUsZGZfelBfVkZFX01MRSxkZl96SV9WRkVfTUxFKQpQX3RydWUgPSBkYXRhLmZyYW1lKHZhbD1Qc2ltLCB0aW1lPXRpbWUpClBfdHJ1ZSRvdXRwdXQgPSAicHJlc3N1cmUgKG1tSGcpIgpJX3RydWUgPSBkYXRhLmZyYW1lKHZhbD1mbG93LCB0aW1lPXRpbWUpCklfdHJ1ZSRvdXRwdXQgPSAiaW5mbG93IChtbC9taW4pIgp0cnVlX291dCA9IHJiaW5kKFBfdHJ1ZSwgSV90cnVlKQoKb2JzUCA9IGRhdGEuZnJhbWUodmFsdWU9ZGF0YV9QSSR5UCwgdGltZSA9IGRhdGFfUEkkdFAsIG91dHB1dCA9ICJwcmVzc3VyZSAobW1IZykiKQpvYnNJID0gZGF0YS5mcmFtZSh2YWx1ZT1kYXRhX1BJJHlJLCB0aW1lID0gZGF0YV9QSSR0SSwgb3V0cHV0ID0gImluZmxvdyAobWwvbWluKSIpCm9icyA9IHJiaW5kKG9ic1Asb2JzSSkKCnBsX3ByZWQ9Z2dwbG90KCkrCiAgZ2VvbV9wb2ludChkYXRhID0gb2JzLCBhZXMoeT12YWx1ZSwgeD10aW1lLCBjb2xvdXIgPSAib2JzZXJ2ZWQiKSwgc2hhcGUgPSAyMCkrCiAgZ2VvbV9saW5lKGRhdGEgPSBwcmVkX2RmLCBhZXMoeT1tZWFuLCB4PXRpbWUsIGxpbmV0eXBlID0gIm1lYW4iKSwgc2l6ZT0wLjkpKwogIGdlb21fbGluZShkYXRhID0gdHJ1ZV9vdXQsIGFlcyh5PXZhbCwgeD10aW1lLCBsaW5ldHlwZT0idHJ1ZSIpLCBzaXplPTAuOSkrCiAgZ2VvbV9yaWJib24oZGF0YSA9IHByZWRfZGYsYWVzKHltaW49bG93ZXIsIHltYXg9dXBwZXIsIHg9dGltZSwgZmlsbCA9ICI5MCUgQ0kiKSwgYWxwaGEgPSAwLjMpKwogIGZhY2V0X2dyaWQob3V0cHV0fm1vZGVsLHNjYWxlcyA9ICJmcmVlIikrCiAgZ2VvbV9wb2ludChkYXRhID0gZGZfeiwgYWVzKHg9eiwgeT15LCBzaGFwZT0iWiIpLCBzaXplPTMpKwogIHNjYWxlX2ZpbGxfbWFudWFsKCIiLHZhbHVlcz1jKCI5MCUgQ0kiID0gImdyZXkxMiIpKSsKICB0aGVtZV9idygpK3hsYWIoInRpbWUgKHNlYykiKSt5bGFiKCIiKSsKICBzY2FsZV9zaGFwZV9tYW51YWwoIiIsIHZhbHVlcyA9IGMoIloiID0gNCkpKwogIHRoZW1lKCNsZWdlbmQucG9zaXRpb24gPSAibm9uZSIsCiAgICAgICAgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpLAogICAgICAgIGxlZ2VuZC5wb3NpdGlvbj0iYm90dG9tIiwKICAgICAgICBzdHJpcC50ZXh0LnggPSBlbGVtZW50X3RleHQoc2l6ZSA9IDEzKSwKICAgICAgICBzdHJpcC50ZXh0LnkgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDEwKSkKKHBsX3ByZWQ9cGxfcHJlZCArIHRoZW1lKHBhbmVsLmdyaWQubWFqb3IgPSBlbGVtZW50X2JsYW5rKCksIHBhbmVsLmdyaWQubWlub3IgPSBlbGVtZW50X2JsYW5rKCkpKQoKZ2dzYXZlKCJmaWd1cmVzL0FwcHJfZGVsdGFfcHJlZF9ub2lzZS5wZGYiLCBwbG90ID0gcGxfcHJlZCwgd2lkdGggPSAxNiwgaGVpZ2h0ID0gMTAsIHVuaXRzID0gImNtIikKYGBgCgoKCmBgYHtyfQpzZXNzaW9uSW5mbygpCmBgYAoKCgo=